Forensic Science MCP US Visualizations
Kyra E Stull, Elaine Y Chu
2022-10-27
- Introduction
- Setup
- Figure 1 – Data Availability
- Figure 2 – Age and sex distributions of the sample separated by the training and testing subsets.
- Figure 5 – Parameter specifications selected for all univariate models and visualized by indicator type.
- Figure 6 – Parameter specifications per variable, separated by indicator type.
- Figure 7 – Univariate Residuals (left) and Absolute Residuals (right), Univariate High Accuracy and Low TMNLP.
- Figure 8 - Dental Development Ordinal CrI
- Figure 9 - Epiphyseal Fusion Ordinal CrI 4-Stage
- Figure 10 - Epiphyseal Fusion Ordinal CI 7-Stage
- Figure 11 - Collapsed and Expanded Comparison
- Figure 12 - Conditional Dependent Long Bone and Dental Model Prediction Comparison
- Figure 13 – Absolute Residuals (top) and Residuals (bottom) Multivariate Conditionally Dependent (Collapsed)
- Figure 14 – Absolute Residuals (top) and Residuals (bottom), Univariate + Conditionally Dependent 18-variable Model
- Figure 15 – KL-Div
- Figure 16 - Comparisons of test predictions for PC_Oss, FBDE_EF, max_M2, and RDL
- Figure 17 - Sex-specific Ordinal CrI Comparison
Introduction
This RMarkdown / HTML file is provided as an exact demonstration of figure creation for the manuscript (journal assigned ID: forensicsci-1982685):
Stull et al. 2023 – Subadult Age Estimation using the Mixed Cumulative Probit and a Contemporary United States Population. Forensic Sciences. Special Issue: Estimating Age in Forensic Anthropology. Special Issue editors: Dr. Kanya Godde and Dr. Rebecca Taylor.
Setup
This chunk of code sets up your R Environment up for success
by:
1. Clearing the global environment of any residual objects
2. Loading pertinent libraries used in the subsequent script
Note: Do not forget to set the folder fs_mcp_us as your working directory!
Within the fs_mcp_us folder should be the following sub-folders that hold all 1) data and 2) resulting files for univariate and multivariate model optimization needed to generate all of the subsequent results and tables provided in the manuscript, and 3) ordinal credible interval tables and model test predictions:
- /ci-tables
- /data
- /results_dent
- /results_ef_oss_c
- /results_eighteen_var
- /results_eighteen_var_c
- /results_lb
- /results_prox_dist
- /results_prox_dist_c
- /results_univariate
- /test-predictions
# Clear global environment
rm(list=ls())
# Load libraries
library(tidyverse)
library(magrittr)
library(ggpubr)
library(yada)
library(doParallel)
# Don't forget to set your working directory!This publication uses a slew of specified figure specifications to standardize figure appearance, color theme, and axis labels, as such:
mytheme <- theme_minimal() +
theme(legend.title=element_text(size=18), legend.text=element_text(size=16),
legend.position="top", axis.title=element_text(size=24),
axis.text=element_text(size=20), strip.text=element_text(size=20),
panel.background=element_rect(fill=NA))
npg_pal <- c("#e64b35ff","#4dbbd5ff","#00a087ff","#3c5488ff",
"#f39b7fff","#8491b4ff","#91d1c2ff","#b09c85ff")
age_lab <- "Age (years)"Figure 1 – Data Availability
Data for Figure 1 and Figure 2 are the reformatted training and
testing samples generated through the MCP Vignette. The
original, full dataset is provided in the data folder as
SVAD_US.csv and the function
caret::createDataPartition was used to separate the data
into a 75% training and 25% testing set with equal proportion of sexes
between the two, and the set.seed() defined as as 719921.
The full code for this process is provided as the script
train-test.R.
Note: The following individuals were removed from the test sample and
figures because they were subsequently identified as typo-outliers after
analyses were run:
SVAD_identifier %in% c(“US_0028”,“US_0390”,“US_0856”,“US_0914”)
train <- read.csv("data/SVAD_US_train_reformatted.csv")
train$model <- "Training"
test <- read.csv("data/SVAD_US_test_reformatted.csv")
test$model <- "Testing"
test %<>% filter(SVAD_identifier != "US_0028", SVAD_identifier != "US_390")
# NOTE: The above individuals were identified as typo-outliers and were
# removed from analyses and visualizations
# Combine training and testing data into a single data frame
data <- rbind(test, train)
data$age_int <- as.integer(data$agey)
data$model <- factor(data$model, levels=c("Training","Testing"))
# Identify unique single-year age cohorts, initialize empty data frame,
# and add column names
age_vec <- sort(unique(data$age_int))
missing_df <- data.frame(matrix(nrow=length(age_vec), ncol=62))
names(missing_df) <- names(data %>% select(FDL:IC_EF))
# Loop through each single-year chronological age and calculate the proportion
# of available data for that age by age indicator, store in missing_df
for(i in age_vec) {
sub <- data %>% filter(age_int==i)
for(j in names(missing_df)) {
N <- nrow(sub[j])
prp <- length(which(is.na(sub[j])))/N
missing_df[i+1,j] <- 1-prp
}
}
# Add single-year chronological age cohorts as a column
missing_df <- cbind(age_int=age_vec, missing_df)# Reshape missing_df for plotting, add new column defining variable type
fig1 <- missing_df %>%
gather(key="var",value="value",-age_int)
fig1$type <- ifelse(grepl("man|max",fig1$var), "dent",
ifelse(grepl("EF|Oss",fig1$var), "ef","lb"))
fig1$var <- factor(fig1$var, levels=names(missing_df))
fig1$type <- factor(fig1$type, levels=c("lb","ef","dent"))
# Plot
ggplot(fig1) + mytheme +
geom_tile(aes(x=age_int, y=var, fill=value*100), color="black") +
theme(legend.position="right", strip.text=element_blank(),
panel.border=element_rect(fill=NA),
panel.spacing.y=unit(0.75,"lines")) +
facet_grid(type~., scales="free") +
labs(x=age_lab, y="Variable", fill="% Available") +
scale_x_continuous(breaks=0:21, expand=c(0,0)) +
scale_y_discrete(limits=rev, expand=c(0,0)) +
scale_fill_gradientn(colors=c("#014f86","#90be6d","#ffea00",
"#f8961e","#f94144"))
Figure 2 – Age and sex distributions of the sample separated by the training and testing subsets.
data %>%
ggplot() + geom_histogram(aes(x=agey, alpha=SEX), position="dodge") +
facet_grid(model ~ .) + scale_alpha_manual(values=c(0.9,0.5)) +
labs(y="Count", x=age_lab, alpha="Sex") + mytheme +
scale_x_continuous(limits=c(0,20))
Figure 5 – Parameter specifications selected for all univariate models and visualized by indicator type.
Data for Figure 5 and Figure 6 use the results of
get_best_univariate_params(), stored in the
results_univariate folder as
“US_all_c_univariate_model_parameters.rds”)
# Import results of get_best_univariate_params()
params <- readRDS("results_univariate/US_all_c_univariate_model_parameters.rds")
# Restructure into new data frame and separate by variable type
var_vec <- unique(params$var)
mod_specs <- params %>% select(var, mean_spec, noise_spec) %>% unique()
mod_specs %<>% mutate(ind_type=ifelse(grepl("man_|max_",var),
"Dental Development",
ifelse(grepl("_EF",var), "Epiphyseal Fusion",
ifelse(grepl("_Oss",var), "Ossification",
"Diaphyseal Dimension"))))
mod_specs$mean_spec <- car::recode(mod_specs$mean_spec,
"'lin_ord'='Linear';
'log_ord'='Logarithmic';
c('pow_law_ord','pow_law')='Power Law'")
mod_specs$noise_spec <- car::recode(mod_specs$noise_spec,
"'const'='Homoskedastic';
'lin_pos_int'='Heteroskedastic'")
mod_specs %<>% mutate(specs=paste(mean_spec, noise_spec, sep=" "))
mod_specs$specs <- factor(mod_specs$specs,
levels=c("Linear Homoskedastic",
"Linear Heteroskedastic",
"Logarithmic Homoskedastic",
"Logarithmic Heteroskedastic",
"Power Law Homoskedastic",
"Power Law Heteroskedastic"))
# Restructure for counts by model specification types
sum_df <- mod_specs %>% group_by(ind_type,mean_spec,noise_spec) %>%
summarize(count=n())
sum_df$ind_type <- factor(sum_df$ind_type, levels=c("Dental Development",
"Epiphyseal Fusion",
"Ossification",
"Diaphyseal Dimension"))
sum_df$noise_spec <- factor(sum_df$noise_spec,
levels=c("Homoskedastic", "Heteroskedastic"))ggplot(sum_df, aes(x=noise_spec, y=mean_spec)) +
geom_point(aes(size=count, color=ind_type),
position=position_dodge2(0.5), pch=19) +
geom_text(aes(label=count, group=ind_type),
position=position_dodge2(0.5), size=10) +
mytheme + scale_size(guide="none") +
labs(x="Noise Specification", y="Mean Specification",
color="Indicator Type") +
scale_color_manual(values=npg_pal[1:4]) +
scale_size(range=c(10,25)) +
guides(color=guide_legend(override.aes = list(size=3)), size="none")
Figure 6 – Parameter specifications per variable, separated by indicator type.
ggplot(mod_specs, aes(x=specs, y=var, color=ind_type)) +
geom_point(pch=19, size=5) + mytheme +
facet_wrap(vars(ind_type), scales="free") +
scale_color_manual(values=npg_pal[1:4]) +
scale_x_discrete(breaks=levels(mod_specs$specs),
labels=c("Linear \n Homoskedastic", "Linear \n Heteroskedastic",
"Logarithmic \n Homoskedastic", "Logarithmic \n Heteroskedastic",
"Power Law \n Homoskedastic", "Power Law \n Heteroskedastic")) +
theme(axis.text.x=element_text(angle=35, hjust=1),
legend.position="none",
strip.background=element_rect(fill="grey85")) +
labs(x="Mean and Noise Specifications", y="Variable")
Figure 7 – Univariate Residuals (left) and Absolute Residuals (right), Univariate High Accuracy and Low TMNLP.
Univariate model test predictions are used in Figure 7, Figure 8,
Figure 15, and Figure 17 to demonstrate differences between variable
type and age estimation using the MCP. These test predictions are housed
in the test-predictions folder and were produced using
univariate_batch_calc().
rdl <- read_csv("test-predictions/RDL_US_all_c_test_predictions.csv")
rdl$model <- "Radius Length"
rdl$resid <- rdl$agey - rdl$xmean
rdl$abs_resid <- abs(rdl$agey - rdl$xmean)
pc_oss <- read_csv("test-predictions/PC_Oss_US_all_c_test_predictions.csv")
pc_oss$model <- "Patella Oss"
pc_oss$resid <- pc_oss$agey - pc_oss$xmean
pc_oss$abs_resid <- abs(pc_oss$agey - pc_oss$xmean)
manP2 <- read_csv("test-predictions/man_PM2_US_all_c_test_predictions.csv")
manP2$model <- "Mandibular PM2"
manP2$resid <- manP2$agey - manP2$xmean
manP2$abs_resid <- abs(manP2$agey - manP2$xmean)
maxM1 <- read_csv("test-predictions/max_M1_US_all_c_test_predictions.csv")
maxM1$model <- "Maxillary M1"
maxM1$resid <- maxM1$agey - maxM1$xmean
maxM1$abs_resid <- abs(maxM1$agey - maxM1$xmean)
fdl <- read_csv("test-predictions/FDL_US_all_c_test_predictions.csv")
fdl$model <- "Femur Length"
fdl$resid <- fdl$agey - fdl$xmean
fdl$abs_resid <- abs(fdl$agey - fdl$xmean)
hpef <- read_csv("test-predictions/HPE_EF_US_all_c_test_predictions.csv")
hpef$model <- "Proximal Humerus EF"
hpef$resid <- hpef$agey - hpef$xmean
hpef$abs_resid <- abs(hpef$agey - hpef$xmean)
tdef <- read_csv("test-predictions/TDE_EF_US_all_c_test_predictions.csv")
tdef$model <- "Distal Tibia EF"
tdef$resid <- tdef$agey - tdef$xmean
tdef$abs_resid <- abs(tdef$agey - tdef$xmean)
uni_models <- rbind(rdl, fdl, pc_oss, manP2, maxM1, hpef, tdef)uni_models %>%
filter(model %in% c("Radius Length","Patella Oss",
"Mandibular PM2", "Maxillary M1", "Proximal Humerus EF", "Distal Tibia EF")) %>%
ggplot(aes(x=agey, y=abs_resid)) +
geom_point(aes(shape=model), color="gray49", alpha=0.75, size=2.8) +
stat_smooth(aes(lty=model, color=model), method="loess", lwd=2, se=F) +
labs(x=age_lab, y="Absolute Residuals (years)", lty="Model",
color="Model", shape="Model") +
mytheme + theme(legend.title=element_blank(),
legend.text=element_text(size=14)) +
scale_y_continuous(limits=c(0,15), breaks=seq(0,15,by=5)) +
scale_shape_manual(values=6:1) +
scale_linetype_manual(values=c(4,1,2,3,5,6)) +
scale_color_manual(values=npg_pal[1:6])
uni_models %>%
filter(model %in% c("Radius Length","Patella Oss",
"Mandibular PM2", "Maxillary M1", "Proximal Humerus EF", "Distal Tibia EF")) %>%
ggplot(aes(x=agey, y=resid)) +
geom_point(aes(shape=model), color="gray49", alpha=0.75, size=2.8) +
stat_smooth(aes(lty=model, color=model), method="loess", lwd=2, se=F) +
labs(x=age_lab, y="Residuals (years)", lty="Model",
color="Model", shape="Model") +
mytheme + theme(legend.title=element_blank(),
legend.text=element_text(size=14)) +
scale_y_continuous(limits=c(-16,16), breaks=seq(-15,15,by=5)) +
scale_shape_manual(values=6:1) +
scale_linetype_manual(values=c(4,1,2,3,5,6)) +
scale_color_manual(values=npg_pal[1:6])
Figure 8 - Dental Development Ordinal CrI
Data for Figure 9 are the combined point and ordinal credible
interval (CrI) data housed in the ci-tables folder and were
generated using generate_ord_ci(). allTeeth is
a dataframe that contains data for all the dental variables.
man_C <- read_csv("ci-tables/ordinal_ci_US_all_c_man_C.csv")
man_I1 <- read_csv("ci-tables/ordinal_ci_US_all_c_man_I1.csv")
man_I2 <- read_csv("ci-tables/ordinal_ci_US_all_c_man_I2.csv")
man_M1 <- read_csv("ci-tables/ordinal_ci_US_all_c_man_M1.csv")
man_M2 <- read_csv("ci-tables/ordinal_ci_US_all_c_man_M2.csv")
man_M3 <- read_csv("ci-tables/ordinal_ci_US_all_c_man_M3.csv")
man_PM1 <- read_csv("ci-tables/ordinal_ci_US_all_c_man_PM1.csv")
man_PM2 <- read_csv("ci-tables/ordinal_ci_US_all_c_man_PM2.csv")
max_C <- read_csv("ci-tables/ordinal_ci_US_all_c_max_C.csv")
max_I1 <- read_csv("ci-tables/ordinal_ci_US_all_c_max_I1.csv")
max_I2 <- read_csv("ci-tables/ordinal_ci_US_all_c_max_I2.csv")
max_M1 <- read_csv("ci-tables/ordinal_ci_US_all_c_max_M1.csv")
max_M2 <- read_csv("ci-tables/ordinal_ci_US_all_c_max_M2.csv")
max_M3 <- read_csv("ci-tables/ordinal_ci_US_all_c_max_M3.csv")
max_PM1 <- read_csv("ci-tables/ordinal_ci_US_all_c_max_PM1.csv")
max_PM2 <- read_csv("ci-tables/ordinal_ci_US_all_c_max_PM2.csv")
man_C$variable <- "Canine"
man_C$jaw <- "Mandibular"
man_I1$variable <- "Central Incisor"
man_I1$jaw <- "Mandibular"
man_I2$variable <- "Lateral Incisor"
man_I2$jaw <- "Mandibular"
man_M1$variable <- "First Molar"
man_M1$jaw <- "Mandibular"
man_M2$variable <- "Second Molar"
man_M2$jaw <- "Mandibular"
man_M3$variable <- "Third Molar"
man_M3$jaw <- "Mandibular"
man_PM1$variable <- "First Premolar"
man_PM1$jaw <- "Mandibular"
man_PM2$variable <- "Second Premolar"
man_PM2$jaw <- "Mandibular"
max_C$variable <- "Canine"
max_C$jaw <- "Maxillary"
max_I1$variable <- "Central Incisor"
max_I1$jaw <- "Maxillary"
max_I2$variable <- "Lateral Incisor"
max_I2$jaw <- "Maxillary"
max_M1$variable <- "First Molar"
max_M1$jaw <- "Maxillary"
max_M2$variable <- "Second Molar"
max_M2$jaw <- "Maxillary"
max_M3$variable <- "Third Molar"
max_M3$jaw <- "Maxillary"
max_PM1$variable <- "First Premolar"
max_PM1$jaw <- "Maxillary"
max_PM2$variable <- "Second Premolar"
max_PM2$jaw <- "Maxillary"
allTeeth <- rbind(man_C, man_I1, man_I2, man_PM1, man_PM2, man_M1, man_M2, man_M3, max_C, max_I1, max_I2, max_C, max_PM1, max_PM2, max_M1, max_M2, max_M3)
allTeeth$variable <- factor(allTeeth$variable)
allTeeth$jaw <- factor(allTeeth$jaw)
allTeeth$ord_stage <- car::recode(allTeeth$ord_stage,
"0='1';1='2';2='3';3='4';4='5';5='6';6='7';
7='8';8='9';9='10';10='11';11='12/13'")
allTeeth$ord_stage <- factor(allTeeth$ord_stage,
levels=c('1','2','3','4','5','6','7','8',
'9','10','11','12/13'))ggplot(allTeeth) +
geom_errorbar(aes(y=variable, xmin=lower95, xmax=upper95, color=jaw),
lwd=2) +
geom_point(aes(x=point_est, y=variable, color=jaw), size=3) +
scale_y_discrete(limits = c("Third Molar","Second Molar", "First Molar",
"Second Premolar", "First Premolar", "Canine",
"Lateral Incisor", "Central Incisor")) +
scale_color_manual(values=c("Mandibular"="black", "Maxillary"="#5f93a1")) +
facet_grid(ord_stage ~ .) + labs(x=age_lab, y="Ordinal Stages per Tooth") +
theme_minimal() + theme(legend.title=element_blank(),
strip.background=element_rect(fill="grey85"),
strip.text=element_text(size=25),
axis.text=element_text(size=25),
axis.title=element_text(size=27),
legend.text=element_text(size=27),
panel.border=element_rect(fill=NA),
panel.spacing.y=unit(0.75,"lines"),
legend.position="top") +
scale_x_continuous(breaks=seq(0,22,by=2))
Figure 9 - Epiphyseal Fusion Ordinal CrI 4-Stage
Data for Figure 10 are also housed in the ci-tables folder
and were generated using generate_ord_ci().
allEF contains the point and 95% CrI for all epiphyseal
fusion variables.
fbde_ef <- read_csv("ci-tables/ordinal_ci_US_all_c_FBDE_EF.csv")
fbpe_ef <- read_csv("ci-tables/ordinal_ci_US_all_c_FBPE_EF.csv")
fde_ef <- read_csv("ci-tables/ordinal_ci_US_all_c_FDE_EF.csv")
fgt_ef <- read_csv("ci-tables/ordinal_ci_US_all_c_FGT_EF.csv")
flt_ef <- read_csv("ci-tables/ordinal_ci_US_all_c_FLT_EF.csv")
fh_ef <- read_csv("ci-tables/ordinal_ci_US_all_c_FH_EF.csv")
hde_ef <- read_csv("ci-tables/ordinal_ci_US_all_c_HDE_EF.csv")
hme_ef <- read_csv("ci-tables/ordinal_ci_US_all_c_HME_EF.csv")
hpe_ef <- read_csv("ci-tables/ordinal_ci_US_all_c_HPE_EF.csv")
ic_ef <- read_csv("ci-tables/ordinal_ci_US_all_c_IC_EF.csv")
ilis_ef <- read_csv("ci-tables/ordinal_ci_US_all_c_ILIS_EF.csv")
ispr_ef <- read_csv("ci-tables/ordinal_ci_US_all_c_ISPR_EF.csv")
rde_ef <- read_csv("ci-tables/ordinal_ci_US_all_c_RPE_EF.csv")
rpe_ef <- read_csv("ci-tables/ordinal_ci_US_all_c_RDE_EF.csv")
ude_ef <- read_csv("ci-tables/ordinal_ci_US_all_c_UDE_EF.csv")
upe_ef <- read_csv("ci-tables/ordinal_ci_US_all_c_UPE_EF.csv")
tde_ef <- read_csv("ci-tables/ordinal_ci_US_all_c_TDE_EF.csv")
tpe_ef <- read_csv("ci-tables/ordinal_ci_US_all_c_TPE_EF.csv")
fbde_ef_7 <- read_csv("ci-tables/ordinal_ci_US_all_FBDE_EF.csv")
fbpe_ef_7 <- read_csv("ci-tables/ordinal_ci_US_all_FBPE_EF.csv")
fde_ef_7 <- read_csv("ci-tables/ordinal_ci_US_all_FDE_EF.csv")
fh_ef_7 <- read_csv("ci-tables/ordinal_ci_US_all_FH_EF.csv")
hde_ef_7 <- read_csv("ci-tables/ordinal_ci_US_all_HDE_EF.csv")
hpe_ef_7 <- read_csv("ci-tables/ordinal_ci_US_all_HPE_EF.csv")
rde_ef_7 <- read_csv("ci-tables/ordinal_ci_US_all_RDE_EF.csv")
rpe_ef_7 <- read_csv("ci-tables/ordinal_ci_US_all_RPE_EF.csv")
tde_ef_7 <- read_csv("ci-tables/ordinal_ci_US_all_TDE_EF.csv")
tpe_ef_7 <- read_csv("ci-tables/ordinal_ci_US_all_TPE_EF.csv")
ude_ef_7 <- read_csv("ci-tables/ordinal_ci_US_all_UDE_EF.csv")
upe_ef_7 <- read_csv("ci-tables/ordinal_ci_US_all_UPE_EF.csv")
hpe_ef_F <- read_csv("ci-tables/ordinal_ci_US_all_c_HPE_EF_sex-F.csv")
hpe_ef_M <- read_csv("ci-tables/ordinal_ci_US_all_c_HPE_EF_sex-M.csv")
tpe_ef_F <- read_csv("ci-tables/ordinal_ci_US_all_c_TPE_EF_sex-F.csv")
tpe_ef_M <- read_csv("ci-tables/ordinal_ci_US_all_c_TPE_EF_sex-M.csv")
fbde_ef$variable <- "Distal Fibula"
fbpe_ef$variable <- "Proximal Fibula"
fde_ef$variable <- "Distal Femur"
fgt_ef$variable <- "Greater Trochanter"
flt_ef$variable <- "Lesser Trochanter"
fh_ef$variable <- "Proximal Femur"
hde_ef$variable <- "Distal Humerus"
hme_ef$variable <- "Medial Epicondyle of the Humerus"
hpe_ef$variable <- "Proximal Humerus"
rpe_ef$variable <- "Proximal Radius"
rde_ef$variable <- "Distal Radius"
upe_ef$variable <- "Proximal Ulna"
ude_ef$variable <- "Distal Ulna"
tpe_ef$variable <- "Proximal Tibia"
tde_ef$variable <- "Distal Tibia"
ilis_ef$variable <- "Ilio-Ischium"
ispr_ef$variable <- "Ischio-pubic Ramus"
ic_ef$variable <- "Iliac Crest"
fbde_ef_7$variable <- "Distal Fibula"
fbpe_ef_7$variable <- "Proximal Fibula"
fde_ef_7$variable <- "Distal Femur"
fh_ef_7$variable <- "Proximal Femur"
hde_ef_7$variable <- "Distal Humerus"
hpe_ef_7$variable <- "Proximal Humerus"
rde_ef_7$variable <- "Distal Radius"
rpe_ef_7$variable <- "Proximal Radius"
tde_ef_7$variable <- "Distal Tibia"
tpe_ef_7$variable <- "Proximal Tibia"
ude_ef_7$variable <- "Distal Ulna"
upe_ef_7$variable <- "Proximal Ulna"
hpe_ef_F$variable <- "Proximal Humerus"
hpe_ef_M$variable <- "Proximal Humerus"
tpe_ef_F$variable <- "Proximal Tibia"
tpe_ef_M$variable <- "Proximal Tibia"
fbde_ef$bone <- "Fibula"
fbpe_ef$bone <- "Fibula"
fde_ef$bone <- "Femur"
fgt_ef$bone <- "Femur"
flt_ef$bone <- "Femur"
fh_ef$bone <- "Femur"
hde_ef$bone <- "Humerus"
hme_ef$bone <- "Humerus"
hpe_ef$bone <- "Humerus"
rpe_ef$bone <- "Radius"
rde_ef$bone <- "Radius"
upe_ef$bone <- "Ulna"
ude_ef$bone <- "Ulna"
tpe_ef$bone <- "Tibia"
tde_ef$bone <- "Tibia"
ilis_ef$bone <- "Os Coxa"
ispr_ef$bone <- "Os Coxa"
ic_ef$bone <- "Os Coxa"
fbde_ef_7$bone <- "Fibula"
fbpe_ef_7$bone <- "Fibula"
fde_ef_7$bone <- "Femur"
fh_ef_7$bone <- "Femur"
hde_ef_7$bone <- "Humerus"
hpe_ef_7$bone <- "Humerus"
rde_ef_7$bone <- "Radius"
rpe_ef_7$bone <- "Radius"
tde_ef_7$bone <- "Tibia"
tpe_ef_7$bone <- "Tibia"
ude_ef_7$bone <- "Ulna"
upe_ef_7$bone <- "Ulna"
hpe_ef_F$bone <- "Humerus"
hpe_ef_M$bone <- "Humerus"
tpe_ef_F$bone <- "Tibia"
tpe_ef_M$bone <- "Tibia"
allEF <- rbind(fbde_ef, fbpe_ef, fde_ef, fgt_ef, flt_ef, fh_ef, hde_ef,
hme_ef, hpe_ef, rpe_ef, rde_ef, upe_ef, ude_ef, tpe_ef,
tde_ef, ilis_ef, ispr_ef, ic_ef)ggplot(allEF %>% filter(bone != "Os Coxa")) +
geom_errorbar(aes(y=reorder(variable, point_est),
xmin=lower95, xmax=upper95, color=bone), lwd=2) +
geom_point(aes(x=point_est, y=variable), size=3) +
facet_grid(ord_stage ~ .) +
labs(x=age_lab, y="Ordinal Stages per Anatomical Site") +
mytheme + theme(legend.position="top", legend.title = element_blank(),
strip.background=element_rect(fill="grey85"),
panel.border=element_rect(fill=NA),
panel.spacing.y=unit(0.75,"lines")) +
scale_x_continuous(breaks=seq(0,22,by=2)) +
scale_color_manual(values=c("black","cyan","gold","deeppink","grey","purple"))
Figure 10 - Epiphyseal Fusion Ordinal CI 7-Stage
Data for Figure 11 are the point and 95% CrI for the proximal and distal epiphyseal fusion variables.
ef_seven <- rbind(fbde_ef_7, fbpe_ef_7, fde_ef_7, fh_ef_7, hde_ef_7, hpe_ef_7,
rpe_ef_7, rde_ef_7, upe_ef_7, ude_ef_7, tpe_ef_7, tde_ef_7)
ef_seven$ord_stage <- car::recode(ef_seven$ord_stage,
"0='0';1='1';2='1/2';3='2';4='2/3';
5='3';6='4'")
ef_seven$ord_stage <- factor(ef_seven$ord_stage,
levels=c('0','1','1/2','2','2/3','3','4'))ggplot(ef_seven) +
geom_errorbar(aes(y=reorder(variable, point_est),
xmin=lower95, xmax=upper95, color=bone), lwd=2) +
geom_point(aes(x=point_est, y=variable), size=3) +
facet_grid(ord_stage ~ .) +
labs(x=age_lab, y="Ordinal Stages per Anatomical Site") +
mytheme + theme(legend.position="top", legend.title = element_blank(),
strip.background=element_rect(fill="grey85"),
panel.border=element_rect(fill=NA),
panel.spacing.y=unit(0.75,"lines")) +
scale_x_continuous(breaks=seq(0,22,by=2)) +
scale_color_manual(values=c("black","cyan","gold","deeppink","grey","purple"))
Figure 11 - Collapsed and Expanded Comparison
Data for Figure 12 contain the collapsed and uncollapsed (expanded) point and 95% CrI for the proximal and distal epiphyseal fusion variables. To make comparisons, the uncollapsed data stages describing fusion (1/2, 2, 2/3) are combined into a mean point estimate and the 95% CrI encompasses the minimum and maximum bounds from the three stages.
c_df <- allEF[c(grep("Proximal",allEF$variable),grep("Distal",allEF$variable)),]
c_df$type <- "Collapsed"
c_df$variable <- car::recode(c_df$variable,"'Distal Fibula'='FBDE_EF';
'Proximal Fibula'='FBPE_EF';
'Distal Femur'='FDE_EF';
'Proximal Femur'='FH_EF';
'Distal Humerus'='HDE_EF';
'Proximal Humerus'='HPE_EF';
'Distal Radius'='RDE_EF';
'Proximal Radius'='RPE_EF';
'Distal Tibia'='TDE_EF';
'Proximal Tibia'='TPE_EF';
'Distal Ulna'='UDE_EF';
'Proximal Ulna'='UPE_EF'")n_df0 <- ef_seven[c(grep("Proximal",ef_seven$variable),
grep("Distal",ef_seven$variable)),]
n_df0$type <- "Uncollapsed"
n_df0$variable <- car::recode(n_df0$variable,"'Distal Fibula'='FBDE_EF';
'Proximal Fibula'='FBPE_EF';
'Distal Femur'='FDE_EF';
'Proximal Femur'='FH_EF';
'Distal Humerus'='HDE_EF';
'Proximal Humerus'='HPE_EF';
'Distal Radius'='RDE_EF';
'Proximal Radius'='RPE_EF';
'Distal Tibia'='TDE_EF';
'Proximal Tibia'='TPE_EF';
'Distal Ulna'='UDE_EF';
'Proximal Ulna'='UPE_EF'")
n_df <- data.frame(ord_stage=c_df$ord_stage)
n_df$point_est <- c(n_df0$point_est[1:2],mean(n_df0$point_est[3:6]),
n_df0$point_est[7:9],mean(n_df0$point_est[10:13]),
n_df0$point_est[14:16],mean(n_df0$point_est[17:20]),
n_df0$point_est[21:23],mean(n_df0$point_est[24:27]),
n_df0$point_est[28:30],mean(n_df0$point_est[31:34]),
n_df0$point_est[35:37],mean(n_df0$point_est[38:41]),
n_df0$point_est[42:44],mean(n_df0$point_est[45:48]),
n_df0$point_est[49:51],mean(n_df0$point_est[52:55]),
n_df0$point_est[56:58],mean(n_df0$point_est[59:62]),
n_df0$point_est[63:65],mean(n_df0$point_est[66:69]),
n_df0$point_est[70:72],mean(n_df0$point_est[73:76]),
n_df0$point_est[77:79],mean(n_df0$point_est[80:83]),
n_df0$point_est[84])
n_df$lower99 <- c(n_df0$lower99[1:2],mean(n_df0$lower99[3:6]),
n_df0$lower99[7:9],mean(n_df0$lower99[10:13]),
n_df0$lower99[14:16],mean(n_df0$lower99[17:20]),
n_df0$lower99[21:23],mean(n_df0$lower99[24:27]),
n_df0$lower99[28:30],mean(n_df0$lower99[31:34]),
n_df0$lower99[35:37],mean(n_df0$lower99[38:41]),
n_df0$lower99[42:44],mean(n_df0$lower99[45:48]),
n_df0$lower99[49:51],mean(n_df0$lower99[52:55]),
n_df0$lower99[56:58],mean(n_df0$lower99[59:62]),
n_df0$lower99[63:65],mean(n_df0$lower99[66:69]),
n_df0$lower99[70:72],mean(n_df0$lower99[73:76]),
n_df0$lower99[77:79],mean(n_df0$lower99[80:83]),
n_df0$lower99[84])
n_df$lower95 <- c(n_df0$lower95[1:2],mean(n_df0$lower95[3:6]),
n_df0$lower95[7:9],mean(n_df0$lower95[10:13]),
n_df0$lower95[14:16],mean(n_df0$lower95[17:20]),
n_df0$lower95[21:23],mean(n_df0$lower95[24:27]),
n_df0$lower95[28:30],mean(n_df0$lower95[31:34]),
n_df0$lower95[35:37],mean(n_df0$lower95[38:41]),
n_df0$lower95[42:44],mean(n_df0$lower95[45:48]),
n_df0$lower95[49:51],mean(n_df0$lower95[52:55]),
n_df0$lower95[56:58],mean(n_df0$lower95[59:62]),
n_df0$lower95[63:65],mean(n_df0$lower95[66:69]),
n_df0$lower95[70:72],mean(n_df0$lower95[73:76]),
n_df0$lower95[77:79],mean(n_df0$lower95[80:83]),
n_df0$lower95[84])
n_df$upper95 <- c(n_df0$upper95[1:2],mean(n_df0$upper95[3:6]),
n_df0$upper95[7:9],mean(n_df0$upper95[10:13]),
n_df0$upper95[14:16],mean(n_df0$upper95[17:20]),
n_df0$upper95[21:23],mean(n_df0$upper95[24:27]),
n_df0$upper95[28:30],mean(n_df0$upper95[31:34]),
n_df0$upper95[35:37],mean(n_df0$upper95[38:41]),
n_df0$upper95[42:44],mean(n_df0$upper95[45:48]),
n_df0$upper95[49:51],mean(n_df0$upper95[52:55]),
n_df0$upper95[56:58],mean(n_df0$upper95[59:62]),
n_df0$upper95[63:65],mean(n_df0$upper95[66:69]),
n_df0$upper95[70:72],mean(n_df0$upper95[73:76]),
n_df0$upper95[77:79],mean(n_df0$upper95[80:83]),
n_df0$upper95[84])
n_df$upper99 <- c(n_df0$upper99[1:2],mean(n_df0$upper99[3:6]),
n_df0$upper99[7:9],mean(n_df0$upper99[10:13]),
n_df0$upper99[14:16],mean(n_df0$upper99[17:20]),
n_df0$upper99[21:23],mean(n_df0$upper99[24:27]),
n_df0$upper99[28:30],mean(n_df0$upper99[31:34]),
n_df0$upper99[35:37],mean(n_df0$upper99[38:41]),
n_df0$upper99[42:44],mean(n_df0$upper99[45:48]),
n_df0$upper99[49:51],mean(n_df0$upper99[52:55]),
n_df0$upper99[56:58],mean(n_df0$upper99[59:62]),
n_df0$upper99[63:65],mean(n_df0$upper99[66:69]),
n_df0$upper99[70:72],mean(n_df0$upper99[73:76]),
n_df0$upper99[77:79],mean(n_df0$upper99[80:83]),
n_df0$upper99[84])
n_df$variable <- c_df$variableggplot(c_df) +
geom_linerange(aes(xmin=lower95, xmax=upper95,
y=ord_stage), color="#4DBBD5FF", size=1) +
geom_point(aes(y=ord_stage, x=lower95, pch="filled", color="blue"), size=3) +
geom_point(aes(y=ord_stage, x=upper95, pch="filled", color="blue"), size=3) +
mytheme + facet_wrap(vars(variable)) +
geom_linerange(data=n_df,
aes(xmin=lower95, xmax=upper95, y=ord_stage),
color="black", lty="dashed", show.legend=F, size=1) +
geom_point(data=n_df,
aes(y=ord_stage, x=lower95, pch="open", color="black"), size=3) +
geom_point(data=n_df,
aes(y=ord_stage, x=upper95, pch="open", color="black"), size=3) +
scale_shape_manual(values=c("filled"=19, "open"=1),
labels=c("4-Stage", "7-Stage")) +
scale_color_manual(values=c("blue"="#4DBBD5FF", "black"="black"),
labels=c("4-Stage", "7-Stage")) +
theme(panel.border=element_rect(fill=NA)) +
labs(y="Ordinal Stage", x=age_lab,
shape=element_blank(), color=element_blank())
Figure 12 - Conditional Dependent Long Bone and Dental Model Prediction Comparison
Data for Figure 13, 14, and 15 include the test predictions from the
conditionally dependent multivariate models. These test predictions are
housed in the test-predictions folder and were generated using
multivariate_batch_calc().
cdep_18UC <- read_csv("test-predictions/cdep_model_US_eighteen_var_test_predictions.csv")
cdep_lbs <- read_csv("test-predictions/cdep_model_US_lb_test_predictions.csv")
cdep_lbs %<>% filter(!(SVAD_identifier %in% c("US_0390","US_0028")))
cdep_ef_ossC <- read_csv("test-predictions/cdep_model_US_ef_oss_c_test_predictions.csv")
cdep_ef_ossC %<>% filter(SVAD_identifier != "US_0856")
cdep_pdC <- read_csv("test-predictions/cdep_model_US_prox_dist_c_test_predictions.csv")
cdep_pdC %<>% filter(SVAD_identifier != "US_0914")
cdep_dent <- read_csv("test-predictions/cdep_model_US_dent_test_predictions.csv")
cdep_18C <- read_csv("test-predictions/cdep_model_US_eighteen_var_c_test_predictions.csv")
cdep_18UC$model <- "C-Dep 18-Var Uncollapsed"
cdep_lbs$model <- "C-Dep Long Bones"
cdep_ef_ossC$model <- "C-Dep EF_Oss"
cdep_pdC$model <- "C-Dep Prox-Dist"
cdep_dent$model <- "C-Dep Dental"
cdep_18C$model <- "C-Dep 18-Var Collapsed"
cdep_18UC$abs_resid <- abs(cdep_18UC$agey - cdep_18UC$xmean)
cdep_18C$abs_resid <- abs(cdep_18C$agey - cdep_18C$xmean)
cdep_lbs$abs_resid <- abs(cdep_lbs$agey - cdep_lbs$xmean)
cdep_dent$abs_resid <- abs(cdep_dent$agey - cdep_dent$xmean)
cdep_pdC$abs_resid <- abs(cdep_pdC$agey - cdep_pdC$xmean)
cdep_ef_ossC$abs_resid <- abs(cdep_ef_ossC$agey - cdep_ef_ossC$xmean)
cdep_18UC$resid <- cdep_18UC$agey - cdep_18UC$xmean
cdep_18C$resid <- cdep_18C$agey - cdep_18C$xmean
cdep_lbs$resid <- cdep_lbs$agey - cdep_lbs$xmean
cdep_pdC$resid <- cdep_pdC$agey - cdep_pdC$xmean
cdep_ef_ossC$resid <- cdep_ef_ossC$agey - cdep_ef_ossC$xmean
cdep_dent$resid <- cdep_dent$agey - cdep_dent$xmean
all_cdep <- rbind(cdep_18UC,cdep_18C,cdep_lbs,cdep_pdC,cdep_ef_ossC,cdep_dent)all_cdep %>% filter(model %in% c("C-Dep Long Bones","C-Dep Dental")) %>%
ggplot(aes(x=agey, color=model)) +
geom_errorbar(aes(ymin=lower95, ymax=upper95, alpha=model), lwd=1.2) +
geom_point(aes(y=xmean), size=1.5) +
geom_abline(slope=1, intercept=0, color="red", linetype="dashed", lwd=1) +
scale_alpha_manual(values=c(1,0.8)) +
labs(x="Known Age (years)", y="Predicted Age (years)",
alpha="Model", color="Model") +
mytheme + scale_color_manual(values=c(npg_pal[6],"black")) +
theme(legend.title=element_blank())
Figure 13 – Absolute Residuals (top) and Residuals (bottom) Multivariate Conditionally Dependent (Collapsed)
all_cdep %>% filter(model != "C-Dep 18-Var Uncollapsed") %>%
ggplot(aes(x=agey, y=abs_resid, shape=model)) +
geom_point(color="gray49", alpha=0.75, size=2.8) +
stat_smooth(aes(lty=model, color=model), se=F, method="loess", lwd=2) +
labs(x=age_lab, y="Absolute Residuals (years)", linetype="Model",
color="Model", shape="Model") + mytheme +
theme(legend.title=element_blank(),
legend.text=element_text(size=14)) +
scale_y_continuous(limits=c(0,8), breaks=seq(0,8,by=2)) +
scale_shape_manual(values=c(1,2,3,4,5),
labels=c("C-Dep 18-Var", "C-Dep Dental",
"C-Dep EF_Oss", "C-Dep Long Bones",
"C-Dep Prox-Dist")) +
scale_linetype_manual(values=c(1,2,3,4,5),
labels=c("C-Dep 18-Var", "C-Dep Dental",
"C-Dep EF_Oss", "C-Dep Long Bones",
"C-Dep Prox-Dist")) +
scale_color_manual(values=c("#0d322c","#1e7662","#a1dceb",
"#ed192e","#28b79f"),
labels=c("C-Dep 18-Var", "C-Dep Dental",
"C-Dep EF_Oss", "C-Dep Long Bones",
"C-Dep Prox-Dist"))
all_cdep %>% filter(model != "C-Dep 18-Var Uncollapsed") %>%
ggplot(aes(x=agey, y=resid, shape=model)) +
geom_point(color="gray49", alpha=0.75, size=2.8) +
stat_smooth(aes(lty=model, color=model), se=F, method="loess", lwd=2) +
labs(x=age_lab, y="Residuals (years)", linetype="Model",
color="Model", shape="Model") + mytheme +
theme(legend.title=element_blank(),
legend.text=element_text(size=14)) +
scale_y_continuous(limits=c(-8,8), breaks=seq(-8,8,by=4)) +
scale_shape_manual(values=c(1,2,3,4,5),
labels=c("C-Dep 18-Var", "C-Dep Dental",
"C-Dep EF_Oss", "C-Dep Long Bones",
"C-Dep Prox-Dist")) +
scale_linetype_manual(values=c(1,2,3,4,5),
labels=c("C-Dep 18-Var", "C-Dep Dental",
"C-Dep EF_Oss", "C-Dep Long Bones",
"C-Dep Prox-Dist")) +
scale_color_manual(values=c("#0d322c","#1e7662","#a1dceb",
"#ed192e","#28b79f"),
labels=c("C-Dep 18-Var", "C-Dep Dental",
"C-Dep EF_Oss", "C-Dep Long Bones",
"C-Dep Prox-Dist"))
Figure 14 – Absolute Residuals (top) and Residuals (bottom), Univariate + Conditionally Dependent 18-variable Model
A few univariate and multivariate models were chosen to compare
between single- and multi-variable model performance. They are combined
into a new dataframe, uni_multi.
maxM2 <- read_csv("test-predictions/max_M2_US_all_c_test_predictions.csv")
maxM2$model <- "Max M2"
maxM2$abs_resid <- abs(maxM2$agey - maxM2$xmean)
maxM2$resid <- maxM2$agey - maxM2$xmean
uni_multi <- rbind(cdep_18C, fdl, maxM2, maxM1, tdef, hpef)uni_multi %>%
ggplot(aes(x=agey, y=abs_resid)) +
geom_point(aes(shape=model), color="gray49", alpha=0.75, size=2.8) +
stat_smooth(aes(lty=model, color=model), method="loess", lwd=2, se=F) +
labs(x=age_lab, y="Absolute Residuals (years)", lty="Model",
color="Model", shape="Model") + mytheme +
theme(legend.title=element_blank(),
legend.text=element_text(size=14)) +
scale_y_continuous(limits=c(0,10), breaks=seq(0,10,by=2)) +
scale_shape_manual(values=c(1,5,2,8,3,4),
labels=c("C-Dep 18-Var","Proximal Humerus EF","Femur Length",
"Max M2","Max M1","Distal Tibia EF")) +
scale_linetype_manual(values=c(1,5,2,6,3,4),
labels=c("C-Dep 18-Var","Proximal Humerus EF","Femur Length",
"Max M2","Max M1","Distal Tibia EF")) +
scale_color_manual(values=c("black","#ed192e","#1e7662",
"#a61120","#28b79f","#a1dceb"),
labels=c("C-Dep 18-Var","Proximal Humerus EF","Femur Length",
"Max M2","Max M1","Distal Tibia EF"))
uni_multi %>%
ggplot(aes(x=agey, y=resid)) +
geom_point(aes(shape=model), color="gray49", alpha=0.75, size=2.8) +
stat_smooth(aes(lty=model, color=model), method="loess", lwd=2, se=F) +
labs(x=age_lab, y="Residuals (years)", lty="Model",
color="Model", shape="Model") + mytheme +
theme(legend.title=element_blank(),
legend.text=element_text(size=14)) +
scale_y_continuous(limits=c(-10,10), breaks=seq(-10,10,by=5)) +
scale_shape_manual(values=c(1,5,2,8,3,4),
labels=c("C-Dep 18-Var","Proximal Humerus EF","Femur Length",
"Max M2","Max M1","Distal Tibia EF")) +
scale_linetype_manual(values=c(1,5,2,6,3,4),
labels=c("C-Dep 18-Var","Proximal Humerus EF","Femur Length",
"Max M2","Max M1","Distal Tibia EF")) +
scale_color_manual(values=c("black","#ed192e","#1e7662",
"#a61120","#28b79f","#a1dceb"),
labels=c("C-Dep 18-Var","Proximal Humerus EF","Femur Length",
"Max M2","Max M1","Distal Tibia EF"))
Figure 15 – KL-Div
Data for Figure 16 contain the Kullback-Leibler Divergence results
calulated using the script kl_div.R.
xcalc <- seq(0,23,by=0.01)
fprior <- readRDS("data/priorx.rds")
US_0088 <- readRDS("data/US_0088.rds")
known_age <- US_0088$ageyNext, because the posterior density distributions and kl-div values were stored in a list, they must be pulled out of said list for use:
kl_list <- readRDS("data/US_88_list.rds")
post_multi_cdep_c <- kl_list$post[[1]]
post_multi_cindep_c <- kl_list$post[[2]]
post_multi_cdep <- kl_list$post[[3]]
post_multi_cindep <- kl_list$post[[4]]
post_lb_cdep <- kl_list$post[[5]]
post_lb_cindep <- kl_list$post[[6]]
post_ef_oss_cdep <- kl_list$post[[7]]
post_ef_oss_cindep <- kl_list$post[[8]]
post_dent_cdep <- kl_list$post[[9]]
post_dent_cindep <- kl_list$post[[10]]
post_rdl_homo <- kl_list$post[[11]]
post_rdl_hetero <- kl_list$post[[12]]
post_fdl_homo <- kl_list$post[[13]]
post_fdl_hetero <- kl_list$post[[14]]
post_man_PM2_homo <- kl_list$post[[15]]
post_man_PM2_hetero <- kl_list$post[[16]]
post_pc_oss_homo <- kl_list$post[[17]]
post_pc_oss_hetero <- kl_list$post[[18]]
kl_div_multi_cdep_c <- kl_list$kl_div[[1]]
kl_div_multi_cindep_c <- kl_list$kl_div[[2]]
kl_div_multi_cdep <- kl_list$kl_div[[3]]
kl_div_multi_cindep <- kl_list$kl_div[[4]]
kl_div_lb_cdep <- kl_list$kl_div[[5]]
kl_div_lb_cindep <- kl_list$kl_div[[6]]
kl_div_ef_oss_cdep <- kl_list$kl_div[[7]]
kl_div_ef_oss_cindep <- kl_list$kl_div[[8]]
kl_div_dent_cdep <- kl_list$kl_div[[9]]
kl_div_dent_cindep <- kl_list$kl_div[[10]]
kl_div_rdl_homo <- kl_list$kl_div[[11]]
kl_div_rdl_hetero <- kl_list$kl_div[[12]]
kl_div_fdl_homo <- kl_list$kl_div[[13]]
kl_div_fdl_hetero <- kl_list$kl_div[[14]]
kl_div_man_PM2_homo <- kl_list$kl_div[[15]]
kl_div_man_PM2_hetero <- kl_list$kl_div[[16]]
kl_div_pc_oss_homo <- kl_list$kl_div[[17]]
kl_div_pc_oss_hetero <- kl_list$kl_div[[18]]data.frame(xcalc=xcalc, multi=post_multi_cdep_c$density,
lb=post_lb_cdep$density,
ef_oss=post_ef_oss_cdep$density,
dent=post_dent_cdep$density,
rdl=post_rdl_homo$density,
pc_oss=post_pc_oss_homo$density) %>%
ggplot() +
geom_line(aes(x=xcalc, y=multi, color="multi"), size=1) +
geom_line(aes(x=xcalc, y=lb, color="lb"), size=1) +
geom_line(aes(x=xcalc, y=ef_oss, color="ef_oss"), size=1) +
geom_line(aes(x=xcalc, y=dent, color="dent"), size=1) +
geom_line(aes(x=xcalc, y=rdl, color="rdl"), lty="dashed", size=1) +
geom_line(aes(x=xcalc, y=pc_oss, color="pc_oss"), lty="dashed", size=1) +
mytheme + geom_vline(xintercept=known_age, color="black", size=1) +
scale_color_manual(values=c("multi"=npg_pal[1],
"lb"=npg_pal[2],
"ef_oss"=npg_pal[3],
"dent"=npg_pal[4],
"rdl"=npg_pal[5],
"pc_oss"=npg_pal[6],
"black"="black"),
labels=c("18-Var / Dep",
"Long Bones / Dep",
"Epiphyseal Fusion / Dep",
"Dental Development / Dep",
"RDL / Homo",
"PC_Oss / Homo",
"Known Age")) +
labs(x=age_lab, y="Density",
color="Model", lty="Model") +
theme(legend.position=c(0.8,0.8),
legend.background=element_rect(),
legend.title=element_blank())
data.frame(xcalc=xcalc,prior=fprior,
cdep=post_lb_cdep$density,
cindep=post_lb_cindep$density,
homo=post_rdl_homo$density,
hetero=post_rdl_hetero$density) %>%
ggplot(aes(x=xcalc)) +
geom_line(aes(y=cdep, color="cdep"), size=1) +
geom_line(aes(y=cindep, color="cindep"), size=1) +
geom_line(aes(y=homo, color="homo"), size=1, lty="dashed") +
geom_line(aes(y=hetero, color="hetero"), size=1, lty="dashed") +
mytheme + geom_vline(xintercept=known_age, color="black", size=1) +
scale_color_manual(values=c("cdep"=npg_pal[1],
"cindep"=npg_pal[2],
"homo"=npg_pal[3],
"hetero"=npg_pal[4],
"black"="black"),
labels=c("LB / C-Dep",
"LB / C-Indep",
"RDL / Homo",
"RDL / Hetero",
"Known Age")) +
labs(x=age_lab, y="Density",
color="Model", lty="Model") +
theme(legend.position=c(0.8,0.8),
legend.background=element_rect(),
legend.title=element_blank())
data.frame(xcalc=xcalc,prior=fprior,
cdep=post_dent_cdep$density,
cindep=post_dent_cindep$density) %>%
ggplot() +
geom_line(aes(x=xcalc, y=prior, color="prior"),
lty="dashed", size=1) +
geom_line(aes(x=xcalc, y=cdep, color="cdep"), size=1) +
geom_line(aes(x=xcalc, y=cindep, color="cindep"), size=1) +
mytheme + geom_vline(xintercept=known_age, color="black", size=1) +
scale_color_manual(values=c("prior"=npg_pal[8],
"cdep"=npg_pal[1],
"cindep"=npg_pal[2],
"black"="black"),
labels=c("Prior",
paste0("C-Dep (KL=",
round(kl_div_dent_cdep,2),")"),
paste0("C-Indep (KL=",
round(kl_div_dent_cindep,2),")"),
"Known Age")) +
scale_linetype_manual(values=c("dashed"),labels=c("Prior")) +
labs(x=age_lab, y="Density (Dental Development)",
color="Model", lty="Model") +
theme(legend.position=c(0.8,0.8),
legend.background=element_rect(),
legend.title=element_blank())
data.frame(xcalc=xcalc,prior=fprior,
cdep=post_multi_cdep_c$density,
cindep=post_multi_cindep_c$density) %>%
ggplot() +
geom_line(aes(x=xcalc, y=prior, color="prior"),
lty="dashed", size=1) +
geom_line(aes(x=xcalc, y=cdep, color="cdep"), size=1) +
geom_line(aes(x=xcalc, y=cindep, color="cindep"), size=1) +
mytheme + geom_vline(xintercept=known_age, size=1) +
scale_color_manual(values=c("prior"=npg_pal[8],
"cdep"=npg_pal[1],
"cindep"=npg_pal[2],
"black"="black"),
labels=c("Prior",
paste0("C-Dep (KL=",
round(kl_div_multi_cdep_c,2),")"),
paste0("C-Indep (KL=",
round(kl_div_multi_cindep_c,2),")"),
"Known Age")) +
scale_linetype_manual(values=c("dashed"),labels=c("Prior")) +
labs(x=age_lab, y="Density (18-Var Collapsed)",
color="Model", lty="Model") +
theme(legend.position=c(0.8,0.8),
legend.background=element_rect(),
legend.title=element_blank())
Figure 16 - Comparisons of test predictions for PC_Oss, FBDE_EF, max_M2, and RDL
ggplot(pc_oss) +
geom_errorbar(aes(x=agey, ymin=lower95, ymax=upper95), color="black") +
geom_point(aes(x=agey, y=xmode)) +
geom_abline(slope=1, intercept=0, color="red", linetype="dashed") +
labs(x="Known Age (years)", y="Predicted Age (years)") +
mytheme
ggplot(tdef) +
geom_errorbar(aes(x=agey, ymin=lower95, ymax=upper95), color="black") +
geom_point(aes(x=agey, y=xmode)) +
geom_abline(slope=1, intercept=0, color="red", linetype="dashed") +
labs(x="Known Age (years)", y="Predicted Age (years)") +
mytheme
ggplot(maxM1) +
geom_errorbar(aes(x=agey, ymin=lower95, ymax=upper95), color="black") +
geom_point(aes(x=agey, y=xmode)) +
geom_abline(slope=1, intercept=0, color="red", linetype="dashed") +
labs(x="Known Age (years)", y="Predicted Age (years)") +
mytheme
ggplot(rdl) +
geom_errorbar(aes(x=agey, ymin=lower95, ymax=upper95), color="black") +
geom_point(aes(x=agey, y=xmode)) +
geom_abline(slope=1, intercept=0, color="red", linetype="dashed") +
labs(x="Known Age (years)", y="Predicted Age (years)") +
mytheme
Figure 17 - Sex-specific Ordinal CrI Comparison
Data for Figure 18 contain the point and 95% CrI for HPE_EF and TPE_EF using pooled, female-only, and male-only datasets.
sex_ord <- rbind(hpe_ef_7 %>% mutate(sex="Pooled"),
hpe_ef_F %>% mutate(sex="Female"),
hpe_ef_M %>% mutate(sex="Male"),
tpe_ef_7 %>% mutate(sex="Pooled"),
tpe_ef_F %>% mutate(sex="Female"),
tpe_ef_M %>% mutate(sex="Male"))
sex_ord$ord_stage <- car::recode(sex_ord$ord_stage,
"0='0';1='1';2='1/2';
3='2';4='2/3';5='3';6='4'")sex_ord %>%
ggplot(aes(y=sex)) +
geom_errorbar(aes(xmin=lower95, xmax=upper95, color=ord_stage), lwd=2) +
geom_point(aes(x=point_est), size=3) + mytheme +
facet_grid(ord_stage ~ variable) +
theme(legend.position="none",
strip.background.y=element_rect(fill="gray85")) +
scale_x_continuous(breaks=c(0,10,20)) +
scale_color_manual(values=npg_pal[1:7]) +
labs(x=age_lab, y="Sample")